Cluster dynamical mean field theory of quantum phases on a honeycomb lattice 
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We have studied the ground state of the half-filled Hubbard model on a honeycomb lattice by 
performing the cluster dynamical mean field theory calculations with exact diagonalization on the 
cluster-impurity solver. Through using elaborate numerical analytic continuation, we identify the 
existence of a 'spin liquid' from the on-site interaction U = Q to Uc (between A.Qt and 4.85f) with 
a smooth crossover correspondingly from the charge fiuctuation dominating phase into the charge 
correlation dominating phase. The semi-metallic state exits only at f/ = 0. We further find that the 
magnetic phase transition at Uc from the 'spin liquid' to the Neel antiferromagnetic Mott insulating 
phase is a first-order quantum phase transition. We also show that the charge fiuctuation plays a 
substantial role on keeping the 'spin liquid' phase against the emergence of a magnetic order. 

PACS numbers: 71.10.-w, 71.27.-(-a, 71.30.+h, 71.10.Fd, 71.10.Pm 



Quantum phase transition is a fascinating physics sub- 
ject, which describes an abrupt change of the ground 
state of a quantum many-body system tuned by a non- 
thermal physical parameter, often accompanied with a 
novel quantum emergence phenomenon As a canoni- 
cal quantum phase transition, the Mott transition, from 
metallic to insulating state tuned by electronic Coulomb 
interaction, is one of the most celebrated and difficult 
problems in condensed matter physics Q . The resultant 
insulating state, namely Mott insulator, usually adopts 
spontaneous symmetry breaking in two and three spa- 
tial dimension to form a long-range antiferromagnetic 
(AFM) order to release the spin entropy due to localized 
electrons. Theoretically, the simplest model to capture 
such physics is the standard one-band half-filled Hub- 
bard model. In the large Coulomb interaction limit this 
model reduces to a standard Heisenberg model with an 
AFM order in its ground state. 

Nevertheless, an insulating ground state without any 
spontaneous symmetry breaking, namely spin liquid, 
may arise if there is frustration d H. Actually the 
spin liquid is a genuine Mott insulating state in a sense 
that it is adiabatically separated from a band insulator. 
Spin liquid has been one of the most intriguing issues in 
condensed matter physics since it was introduced nearly 
forty years ago and continuously in intense research 
since it was further proposed to be a parent phase to 
likely lead to high Tc superconductivity However, a 
spin liquid had not been verified for a two-dimensional 
(2D) standard Hubbard or Heisenberg model until a re- 
cent quantum Monte Carlo simulation was done for cor- 
related fermions on a honeycomb lattice f(|] . Through fi- 
nite size extrapolation the simulation surprisingly shows 
that a spin liquid emerges between the semi-metallic and 
AFM Mott insulating states. This has stimulated a new 
surge of strong discussion on the nature of such a quan- 
tum phase and the related phase diagram 0, Q • 

To help to clarify this issue, we performed the cluster 
dynamical mean field theory calculations for the half- 
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FIG. 1: (Color online) Schematic phase diagram of the one- 
band Hubbard model on a honeycomb lattice at half-filling. 
Uc is between A.Qt and 4.85t. 



filled Hubbard model defined on a honeycomb lattice as. 
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where {cia) is the electron creation (annihilation) op- 
erator with spin a or \,) at lattice site i, (ij) represents 
the summation over the nearest neighbors, i > is the 
nearest neighbor hopping integral, and n^o- = cj^Ci^ with 
the on-site Coulomb repulsion U. 

Our calculated results are schematically summarized 
in Fig. [H By the calculations, we show that a disordered 
phase of 'spin liquid' exists from U = to Uc, at which 
it transforms into the Neel AFM insulating phase via a 
first-order quantum phase transition. 

The dynamical mean field theory (DMFT) maps a 
quantum lattice model onto a single lattice site, a quan- 
tum impurity, dynamically coupled to a self-consistently 
determined bath of free electrons that represents the 
rest of the lattice [gj. Thus the DMFT fully consid- 
ers local quantum dynamical fiuctuations, much beyond 
conventional mean-field methods. The DMFT has sub- 
stantially improved our understanding on the nonper- 
turbative properties of correlated electron systems, par- 
ticularly the Mott transition. The cluster dynamical 
mean field theory (CDMFT) is a natural extension of the 
DMFT to include the missed short-ranged spatial corre- 
lations through a proper replacement of a single-site im- 
purity by a cluster of lattice sites, which is constructed 
to reflect the lattice symmetry and local lattice structure 
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features \mUM- The CDMFT has been successfully ap- 
plied to study a variety of ordered phases, and opens an 
avenue to directly study quantum phase transitions. 

Here it should be addressed that for the (C)DMFT the 
thermodynamic limit is naturally taken from the outset 
through a self-consistent procedure 0, . As a nonper- 
turbative approach to treat many-body correlation ef- 
fects, the (C)DMFT works well in the whole coupling 
regime and becomes exact in the two contrary limits of 
both noninteracting and infinite-interacting cases. For 
finite-size quantum Monte Carlo simulations or exact 
diagonalizations, in contrast, the thermodynamic limit 
is extrapolated through finite-size scaling. Correspond- 
ingly, the CDMFT allows for spontaneous symmetry 
breaking, while the finite-size approaches have difhculty 
in finding a long-range order and related phase transition 
or underestimate ordered phases. Thus the CDMFT and 
finite-size approaches are complementary to each other, 
both of which together can give more conclusive results 
than they alone. 

In CDMFT calculations, the target in the self- 
consistent procedure is to obtain a lattice imaginary 
frequency local Green's function matrix Gij{iuJn) (sub- 
scripts i and j being site indices of a chosen cluster) by as- 
suming its self-energy matrix identified as the one of the 
corresponding cluster-impurity Green's function matrix 
GlJ^^ (iujn), derived from the Dyson equation. In order 
to study a ground state, we apply exact diagonalization 
rather than quantum Monte Carlo simulation to solve 
a cluster- impurity model [13, Specifically, we em- 

ted the robust Krylov-Schur algorithm based SLEPc 
to accomplish the large-scale sparse matrix diagonal- 
ization efficiently and stably Here we particularly 
emphasize that to carry out an elaborate numerical ana- 
lytic continuation from an imaginary frequency Green's 
function Gii(iw„) onto a real frequency retarded Green's 
function Gu{uj -\- iO"*") is crucial to unambiguously iden- 
tify whether or not an energy gap exists at a small on-site 
interaction U by checking the density of states 
(DOS) equal to —^ImGii{uj + iO+). In such a way, the 
energy gap resolution can be achieved as high as 10~^i, 
which is hardly reached by other methods. 

We first performed the CDMFT calculations for one- 
dimensional (ID) Hubbard model at half-filling with a 
cluster-impurity model respectively containing two and 
four lattice impurity sites, which serves as a benchmark 
for the further calculations. The calculated results are 
reported in Fig. [21 in quantitative comparison with the 
Bethe Ansatz exact solution [i^. As we see, the two- 
site CDMFT result has been already in excellent agree- 
ment with the exact solution. Especially, by using the 
numerical analytic continuation, we can unambiguously 
identify a finite energy gap immediately develops once 
U is nonzero. Thus the short-ranged spatial correla- 
tions play a dominant role in local spectral functions and 
properties even though the spatial correlations have long- 




FIG. 2: (Color online) Energy gap A, namely single-particle 
spectral gap, as a function of the on-site interaction U in the 
one-dimensional half-filled Hubbard model, calculated with 
the two-site [Nc=2) and four-site {Nc=A) clusters, respec- 
tively. The dashed curve denotes the exact one from the Bethe 
ansatz solution. The inset zooms in the small-C/ dependence. 
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FIG. 3: (Color online) (a) Cluster-impurity model configura- 
tion for a honeycomb lattice. Filled circles denote the impu- 
rity sites. Unfilled circles denote the bath levels. Links rep- 
resent the hopping paths, (b) Density of states of the model 
when U = Q. The blue curve is the exact one. The red one 
was obtained by the numerical analytic continuation. 



ranged power-law behavior in ID Hubbard model |21|. 
It is also well-known that the quantum fiuctuations are 
much stronger in one-dimension than in higher dimen- 
sions. Hence it is highly expected that the higher dimen- 
sional CDMFT results are more reliable and encouraging 
than the one-dimensional ones. 

Figure ISja) schematically shows the cluster-impurity 
model constructed for a honeycomb lattice. Such a model 
reflects sixfold rotational symmetry and an impurity site 
with a lattice coordination number of 3, which are the 
essential features of a honeycomb lattice. We also add a 
direct link between each pair of the nearest bath levels so 
that we can include the propagation of an electron from 
one impurity site through the outside of the cluster (the 
bath) to any other impurity site in the calculations. 

When the on-site interaction U = 0, the half-filled 
Hubbard model on a honeycomb lattice reduces to a set of 
free Dirac fermions with a linear DOS around the Fermi 
energy. As shown in Fig. EJb), the numerical analytic 
continuation can well repeat the exact DOS. Particularly, 
around the Fermi energy both are the strictly same even 
though the DOS is not differentiable at the Fermi energy. 
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FIG. 4: (Color online) Calculated magnetization as a function 
of the on-site Coulomb interaction i7 on a honeycomb lattice. 




FIG. 5: (Color online) Calculated density of states for several 
small values of the on-site Coulomb interaction [/ on a hon- 
eycomb lattice, zooming in around the Fermi energy a; = 0. 
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FIG. 6: (Color online) Calculated density of states for [7 = 4/ 
(red lines) and (7 = 6/ (blue and green lines) on a honeycomb 
lattice, respectively. The dotted curve is the exact DOS at 
f/ = for a reference and the Fermi energy sets to zero. 




In this study, we define the magnetization m=((n^ — 
ni)/2) on a site. Being bipartite, a honeycomb lattice 
can be divided into two sublattices A and B. If a Neel 
AFM state appears, m will alternatively take positive 
and negative along with sublattices A and B, namely 
being staggered magnetization. As we see from Fig. SI 
when U < Uc2 — 4.85/:, a paramagnetic solution of m = 
is stable, namely no spin polarization on each site, but 
over Uc2 this solution is no longer stable and \m\ abruptly 
jumps over 0.2. On the other hand, when U > Uci — 4.6t 
a staggered magnetization solution with |m| > 0.16 is 
stable, namely a Neel AFM phase takes over, but below 
Uci m immediately plummets to zero. Between Ud and 
Uc2, these two solutions coexist. Such a hysteresis be- 
havior indicates this magnetic transition is a first-order 
quantum phase transition. 

We calculated the DOS for a small on-site interaction 
U with extreme caution through elaborate numerical an- 
alytic continuation [ll|. The calculated DOS are then 
plotted in Fig. [5j Similar to the case of ID Hubbard 
model, what we find is that there is also a definite en- 
ergy gap opening at the Fermi energy once U is nonzero. 
In comparison with the case of C/ = by checking the 
enclosed area, it is further shown that the corresponding 
states nearby the Fermi energy are clearly moved away 
from the gap rather than pushed to the two sides of the 
gap. To be specific, for U = OAt, 0.8t, and 1.2/;, the en- 
ergy gap is found as small as A = 2.5 x 10^^/, 1.0 x 10^^/, 
and 0.023/, respectively. For a rather large U, a rela- 
tively large energy gap opens with a substantial portion 



FIG. 7: (Color online) Calculated energy gap A as a function 
of the on-site Coulomb interaction (7 on a honeycomb lattice. 
The inset zooms in the small-f/ dependence. Note that the 
energy gap is in linear i7-dependence after the transition. 



of states moved away from the Fermi energy into below 
—3/ and above 3t, corresponding to the Hubbard band 
states. For a U further larger than Uc2, the system trans- 
forms into the Neel AFM phase. The on-site spin degen- 
eracy is then lifted. The spin-up (spin-down) resolved 
DOS at A sublattice is the same as the spin-down (spin- 
up) resolved DOS at B sublattice, as shown in Fig. |6l 

Figure [7] shows the energy gap A as a function of 
the on-site interaction U, extracted from the calculated 
DOS. For U < 1.6/, the energy gap increases very 
slowly with U, and the function can be represented as 
A// — 0.015(C///)^. After 1.6/, the energy gap increasing 
becomes fast with U increasing. Between Ud and Uc2, 
the [/-dependence of the energy gap shows a hysteresis 
behavior with a sudden change of ^^0.45/, correspond- 
ing to the first-order quantum phase transition. Thus a 
nonzero U definitely induces an energy gap and makes 
the system be in an insulating phase. 

We also calculated the equal time spin-spin correla- 
tion functions (sfs|) and occupancy-occupancy correla- 
tion functions (nirij), with i, j being a pair of the nearest 
sites, next nearest sites, and next next nearest sites, re- 
spectively. We find that all these correlation functions 
keep the sixfold rotational symmetry as the honeycomb 
lattice does before the magnetic transition. This excludes 
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FIG. 8: (Color online) Calculated double occupancy D as a 
function of the on-site Coulomb interaction U ona honeycomb 
lattice. The inset zooms in the hysteresis loop. 

spin valence bond structures, spin anisotropic structures, 
and charge density waves. We can also exclude another 
possibility to break the translation symmetry that ev- 
ery six sites form a benzene-ring-like plaquette to further 
form a plaquette-singiet valence bond solid pattern in a 
honeycomb lattice. If such a pattern exists, the states 
nearby the Fermi energy will be pushed to the two sides 
of the gap rather than moved away. Thus this nonmag- 
netic insulating phase from U = to Uc (between A.6t 
and A.85t) can be classified as a 'spin liquid' phase in the 
sense that it is tuned on by the on-site interaction U. 

To understand the underlying physics, we examine the 
double occupancy, defined as D = (n-^n^) on a site. The 
ground state energy per site Eg = {H)/N of Hamilto- 
nian ([1} is a function of the on-site interaction U . Its 
derivative is nothing but the double occupancy, namely 
dEg/dU = (rifni). Thus the double occupancy D di- 
rectly describes a quantum phase transition tuned by U. 
In Fig. |51 the U -dependence of the D likewise shows a 
hysteresis behavior between Ud and Uc2- This means 
the energy level crossing in the ground state through the 
magnetic transition as U increasing, being a characteris- 
tic of a first-order quantum phase transition [1]. 

The on-site interaction U tunes or controls the Hamil- 
tonian (1) through the double occupancy D. Moreover, 
the localization degree of an electron, as well as the local 
correlation effect, can be quantitatively described by the 
double occupancy. At half-filling, D is between 0.25 and 
0, corresponding respectively to full delocalization and 
complete localization. In addition, the magnetic moment 
is obtained by (ml) = {{2S^f) = 1-2D. 

It is commonly thought that the Mott transition is 
driven by the strong local correlation effect due to the 
on-site interaction U, which is marked by a vanishing or 
very small double occupancy with a large local moment 
on a site 0, [l^]. In contrast, for the Hubbard model on 
a honeycomb lattice, a small U can immediately induce 
a small energy gap opening to tune the system into an 
insulating phase with a large double occupancy namely 
large charge fluctuation, as shown in Fig. [5] The cal- 
culations show that the small- [/ induced energy gap is 



a consequence of the interplay between the zero-DOS at 
the Fermi energy (Dirac Cone band) and local charge 
correlation, not a conventional correlation-driven Mott 
insulating gap. On the other hand, the correlation ef- 
fect will become dominating nearby the magnetic tran- 
sition. Actually Fig. [7] has shown that the energy gap 
becomes a linear function of U after the transition, which 
is the canonical behavior of a correlation-driven Mott in- 
sulator. Thus the 'spin liquid' states nearby U = are 
of charge fluctuation dominating while those nearby Uc 
are of charge correlation dominating, and corresponding 
to the ones found by the quantum Monte Carlo simula- 
tions reported in Ref. @. Nevertheless our calculations 
show that a smooth crossover connects these two contrary 
parts. Meanwhile it is also indicated that the charge 
fluctuation plays a substantial role on keeping the 'spin 
liquid' phase against the emergence of an AFM order. 

In summary, we have performed the cluster dynamical 
mean field theory calculations, allowing for the sponta- 
neous symmetry breaking, to study the ground state of 
the half-filled Hubbard model on a honeycomb lattice. 
We find that a 'spin liquid' exists from U = to about 
A.6t, in which the system takes a smooth crossover corre- 
spondingly from the charge fluctuation dominating phase 
into the charge correlation dominating phase, then it fur- 
ther transforms into the Neel AFM Mott insulating phase 
via a first-order quantum phase transition. 
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